Prom vortex molecules to the Abrikosov lattice in thin mesoscopic superconducting 

disks 
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Stable vortex states are studied in large superconducting thin disks (for numerical purposes we 
considered with radius R — 50£). Configurations containing more than 700 vortices were obtained 
using two different approaches: the nonlinear Ginzburg-Landau (GL) theory and the London ap- 
proximation. To obtain better agreement with results from the GL theory we generalized the London 
theory by including the spatial variation of the order parameter following Clem's ansatz. We find 
that configurations calculated in the London limit are also stable within the Ginzburg-Landau the- 
ory for up to ~ 230 vortices. For large values of the vorticity (typically, L > 100), the vortices 
are arranged in an Abrikosov lattice in the center of the disk, which is surrounded by at least two 
circular shells of vortices. A Voronoi construction is used to identify the defects present in the 
ground state vortex configurations. Such defects cluster near the edge of the disk, but for large L 
also grain boundaries are found which extend up to the center of the disk. 

PACS numbers: 74.20.De, 74.25. Dw, 74.25. Ha 



I. INTRODUCTION 

Vortices appear in several branches of physics, such 
as fluid dynamics^ superfluidity, 2 Bose-Einstein (B-E) 
condensates ) 3 ' 4 ! 5 and superconductivity^! The vortex is 
usually described by a field (for instance, the velocity 
field) which diverges as r -1 as one approaches its coie& 
They can be treated as quasiparticles, since they can be 
created or destroyed, they interact with each other and 
with the interfaces. Unlike in fluid dynamics, in super- 
fluids (including here superconductors and B-E conden- 
sates) vortices are quantized objects. In superconductors, 
for example, they carry a magnetic flux which is a multi- 
ple of the flux quantum, <I>o = hc/2e, and are character- 
ized by a core of area £ 2 - where the superconductivity 
is highly depreciated - surrounded by superconducting 
currents (screened at distances of order A) . Here, £ is the 
coherence length. They have been intensively studied, 
since Abrikosov 6 predicted their existence from the solu- 
tion of the Ginzburg-Landau (GL) equations in a type-II 
superconductor for H& < H < H C 2 (see also Refs. |2l 
and llfib In an infinite, and defect free superconductor, 
vortices arrange themselves in an hexagonal (Abrikosov) 
lattice. 

A detailed phenomenological description of the super- 
conducting state can be derived from the GL theoryjil by 
means of two parameters: the complex order parameter, 
\E', which is related to the superconducting electron den- 
sity, and the vector potential, A. For 7J cl < H <C H c2 , 
each vortex can be viewed as a particle, since inter- vortex 
separations, a, are such that £ <C a ~ A - assuring that 
vortex cores do not overlap - and the major role between 
vortex- vortex interactions is played by the superconduct- 
ing shielding currents. In such cases the London limit 
turns out to be a good approximation of the GL theory, 
beco ming better for higher values of k (see for example 
Refs. l7ll2ll3ll4F) . In this approximation, the supercon- 



ducting electron density is considered constant through- 
out the entire superconductor and the vortex cores are 
represented by singularities in the phase of the order pa- 
rameter. This allows to treat vortices as particles. 

In a thin film of thickness d, the effective magnetic 
field shielding length turns out to be the effective pen- 
etration depth, A = A 2 /d, instead of A~ At distances 
r <C A the electromagnetic interaction is still logarith- 
mic, as in the three dimensional case, but with screening 
length A [However the perpendicular magnetic field and 
the shielding currents decay as r -3 and r~ 2 far away from 
the vortex core for r -C A, instead as exp(— r/A) in the 
bulk case.] Similarly as in the bulk case, in a thin film 
vortices also form an hexagonal Abrikosov latticed 

In mesoscopic superconductors both the geometry and 
size of the specimen influence the vortex configurations, 
due to the interaction between vortices and the surface. 
Therefore, for small enough samples (with sizes compara- 
ble to £), the conventional hexagonal lattice predicted by 
Abrikosov no longer exists, and vortex configurations ad- 
just to the sample geometry, yielding some kind of vortex 
molecule statesiiiiiiiiS* 1 ^ 2 ^ For example, vortices arrange 
themselves in ringlike structures in disks with radii (R) a 
few times £,18,19,20,21,22,23,24,25,26,27,28,29,30 Sudl patterns 

show similarities to what is observed in electrons in arti- 
ficial atoms, where particles obey specific rules for shell 
filling and exhibit magic numbers. Nevertheless when 
overlapping of vortices starts to take place, discrepan- 
cies between vortices and a picture based on particles 
arise, such as the formation of giant vortex states. Also, 
vortex-antivortex configurations may become possible for 
non circular geometries 

Within the London limit the vortex interaction poten- 
tial in a thin disk of arbitrary radius was first calculated 
by Fetteo 3 *. Also in the London limit, vortex configura- 
tions up to L — 8 were studied by Buzdin and Brisor. 35 
for A ^> R (where demagnetization effects can be ne- 
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glected). In the latter limit it is possible to substitute 
the interaction between the vortices and the disk border 
by the interaction between vortices and their images (see 
also Ref . l36j) . Within the London limit one is able to find 
analytical expressions for the energy and forces of an ar- 
bitrary arrangement of vortices inside the disk, since vor- 
tices can be treated as particles. Vortices considered as 
particles were also studied by Monte Carlo and Molecular 
Dynamics simulations. In Ref. |33 vortex configurations 
with up to 2000 vortices were studied and an hexagonal 
lattice was found for thin disks, although they did not 
consider the vortex interaction with the disk edge. Vor- 
tex molecules in long cylinders with radius much larger 
than A were studied by Venegas and Sardella. 38 Other 
geometries were investigated in Refsi^ii^ ., for example. 

In this paper we will study multivortex states where 
many vortices nucleate, yielding a triangular lattice in 
the center of the disk and a ringlike structure close 
to the edges. Within the GL framework several other 
works have been reported regarding vortex states in thin 
disks, 18 i 19 i 21 . 22 i 23 i 24 i 25 i 26 i 27 i 28 i 29 i 30 but they were limited 
to much smaller disk radius. In such small systems the 
formation of multivortex states with high vorticity is not 
allowed and, consequently, it was not possible to study 
the transition from a ringlike structure to an Abrikosov 
lattice, which is the subject of the present paper. 

This paper is organized as follows. The theoretical 
approach is described in Section [H] In Section IIIII low 
vorticity states obtained within the GL and the London 
frameworks are compared. In Sections IIVI and Ivl config- 
urations with up to 700 vortices are investigated, respec- 
tively, by showing the existence of an Abrikosov lattice 
in the center of the disk and by examining the role of 
topological defects in the lattice in order to adjust the 
hexagonal lattice to the radial symmetry close to the disk 
edge. Surface superconductivity in the R — 50£ disk is 
briefly analyzed in Section lVfl Our conclusions are given 
in Section IVTT1 



II. THEORETICAL APPROACH 

For our numerical calculation we used a thin disk of ra- 
dius R — 50£ and thickness d, in which A = \ 2 /d^> 
d, surrounded by vacuum and in the presence of a uni- 
form perpendicular magnetic field Hq. In this regime, 
the demagnetization effects can be neglected, allowing 
one to assume H w Hq. Vortex states in mesoscopic thin 
disks were investigated by us using both the Ginzburg- 
Landau (GL) theory and the London approximation with 
the London gauge V • A = 0. Dimensionless variables are 
used, i.e., the distance is measured in units of the co- 
herence length £, the vector potential in cfi./2e£ and the 
magnetic field in H C 2 — cft/2e£ 2 = k^/2H c . The average 
energy density is written in units of H 2 /8tt (we shall refer 
to it as simply the energy of the system). Also, the vortic- 
ity or the number of vortices in the system will be denoted 
by L (an analogue to the total angular momentum) ,22^ 



Moreover, whenever the distinction among different con- 
figurations with the same L would be necessary, we use 
the notation presented in Ref. |2£j to denote the vortex 
configurations, e.g., for L = 6, (1,5) means 1 vortex in 
the center with 5 around it, and (6) represents 6 vortices 
with none of them in the center of the disk. 

In the framework of the GL theory, the GL equa- 
tions are solved numerically according to the approach 
of Schweigert and Peetersi 2 ^ 2 ^ As we are in the limit 
(d <C £,A), the Ginzburg-Landau equations can be av- 
eraged over the disk thickness, leading to the following 
system of equations, 
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where the supercurrent density is defined by the follow- 



ing, 
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S(z)\^\ 2 (V 2 dO-A)=S(z)\^\ 2 U. (3) 



Above, the superconducting wavefunction, * = |\l/|e 4e , 
satisfies the boundary conditions (— iV 2 D — A) = 
normal to the sample surface and A = A = ^H pcj) 

(since demagnetization effects can be neglected). Here <p 
is the unit vector in the azimuthal direction. The indices 
2D, 3D refer to two- and three-dimensional operators, re- 
spectively. The dimensionless GL energy density is given 
by 



G = Gc 



(4a) 



where 
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-2 |*f + |*| 4 + 2 (V 2D |*|) 2 1 dV, (4b) 



2|*| 2 n 2 + 2k 2 (H - HqY dV, (4c) 



are the core and the electromagnetic energies, respec- 
tively, and the integrations are to be performed over the 
sample volume V. As demagnetization effects can be dis- 
regarded, the above equation reduces to 



(5) 



which was actually the expression used to compute the 
energy of the vortex configurations within the GL the- 
ory. For now on the symbol V will be used for the two- 
dimensional gradient operator. 

The system of Eqs. <|1I2[) were solved by using the ap- 
proach of Ref. 24 for circular disks. A finite-difference 
representation for the order parameter is used on an uni- 
form 2D square grid (x,y), with typically 512 x 512 grid 
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points for the area of the superconductor, which allows 
to have at least 5 grid points inside a length of the order 
of £. We also use the link variable approach^ and an it- 
eration procedure based on the Gauss-Seidel technique to 
find , f. Starting from different randomly generated ini- 
tial conditions and at some specified magnetic field, the 
steady-state solutions of Eqs. (11121) yield different vortex 
configurations, either stable or meta-stable states. 

In the London approximation, the order parameter is 
considered uniform throughout the disk, except for small 
regions with areas of the order of £ 2 , where it drops to 
zero. This can only be accomplished when k 1. Then 
the energy of the system is purely electromagnetic and it 
is given by the sum of the supercurrent and the magnetic 
field energies 



Gl 



2k 2 r 
v J 



dV 



(H-H ) 2 



(6) 



Notice that this expression is a particular case of Eq. ijlc|) 
which is obtained by putting l^l 2 = I everywhere inside 
the disk. In the presence of L vortices, situated at pi 
{i = I, 2, L}, the London equation can be written as 



J 



where 



(7) 



(8) 



i=i 



with pi — (xi, yi) the position of the vortices, J = 
Jq dzj ~ jrf, and = 4>/p- The vortex images at 

(R/ Pi) 2 Pi appear in Eq. (|SJ in order to fulfill the bound- 
ary condition 35 J(R)-p = 0. Instead of writing Eq. Q for 
the vector «7, one may use the streamline function, g(p), 
related to the supercurrent by J = V x (zg) (g(p) can 
be regarded as a local magnetization in the thin film. 42 ) 
At the boundary g(R, <fi) = const, but, as the value of 
this constant is arbitrary, one can impose g(R, 4>) = 0. 
Therefore, Eqs. J7J) and J8J can be expressed as, 
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{R/ Pj ) 2 Pj \ Pj 



\P-Pj\ 



R 



(9) 



Notice that Eq. (JJJ) can also be understood as the limit- 
ing case of the GL equations if one considers \*S?\ = 1 and 
V# = v . Therefore, while vortices are well apart from 
each other (and also the boundary), there exists a rela- 
tion between the streamline function defined above and 
the phase of the order parameter in the GL theory, i.e., 
one can define a complex function of which the real and 
imaginary parts are proportional to g(p) and 

Since in our case (A = A 2 /d> ( > d), demagnetiza- 
tion effects can be neglected 2 ^ and one may write Eq. © 



as 
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2-K^gipi) - H J d 2 pg(p) 



(10) 



where the integration is performed along the thin film 
plane, z = 0. Substituting Eq. JUJl in this formula, and 
after some algebraic manipulation, the London energy is 
expressed by 
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(11) 



where we used = pi/R to simplify the notation. 

The divergence in Eq. (|II|) can be removed by con- 
sidering a cut-off, in which for i = j — > \p+ — pj\ = a£ 
(in not normalized units) and a is a constant. The final 
expression for the London energy can be written as 



where 




self 



^shield 



field 



clf =(|) lnfl-/7 



,(12a) 



(12b) 



is the interaction energy between the i th vortex and the 
radial boundary of the superconductor, 



^shield 



-2H (1 - r? 



(12c) 



represents the interaction between the z th vortex and the 
shielding currents, and 



y = ' it hl 



(rirj) 2 - 2r t ■ r,- 



1 
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(12d) 



is the repulsive energy between vortices i and j. Finally, 



(2/R) 2 Nln{R/a) and 



Jeld 



— R Hn are the en- 



ergies associated with the vortex cores and the external 
magnetic field, respectively. 

Notice that G l allows one to treat the vortices as par- 
ticles. Therefore, simulation techniques appropriate for 
systems of particles may be performed in order to find, 
for example, the ground state of the systemi^ii 4 ^ 5 . In this 
sense, the vortex system behaves (in the London approx- 
imation) similar to a two dimensional system composed 
of equally charged particles interacting through a repul- 
sive logarithmic potential placed in parabolic potential 
wellJiiii Nevertheless, there is a fundamental difference 
between these two systems: the vortex system is confined 
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to the disk of radius R and the influence of the surface 
on the energy is clear from the terms containing vortex 
images, i.e., e^ olf and e^. Notice also that e coro arises 
from the cut-off procedure and is therefore strongly de- 
pendent on the cut-off value at; (we adopted a = 1 in the 
results shown below). The actual energy associated with 
vortex cores and with the spatial variation of the super- 
conducting electron density (\ip(p)\ 2 ) should be evaluated 
by using the GL theory. 

A thin disk with L vortices was simulated by using 
Eq. I12al) . To investigate (meta-)stable states close to the 
equilibrium, we employed a procedure similar to the one 
described in Ref. |4J. First L' vortices were distributed 
randomly inside the disk. Then, a Monte Carlo (MC) 
technique was used to make the system wander in the 
configurational space and arrive at a neighborhood of 
some minimum of Gl- After typically 10 4 MC steps, we 
perform a molecular dynamics (MD) simulation starting 
from the final MC configuration. The final (meta)-stable 
state is achieved after about 10 6 MD steps. In order 
to find the ground state (or states with energies very 
close to it) this trial procedure was repeated more than 
1000 times, each time starting with a different random 
distribution of L' vortices at a given magnetic field Hq. 

To implement the MD we time integrated the Bardeen- 
Stephen equation of motion 4 ^ 



dpi 
alt 



(13) 



where i is the label of the i vortex, 77 is the viscous 
drag coefficient 77 ~ &oH C 2/ p n c 2 (where p n is the nor- 
mal state resistivity). The forces acting on each vortex 
were obtained from — V 'kGh{pi, Pj), where Gl is given by 
Eq. (|12afl and — X7 k is the gradient with respect to the 
coordinate pk- This yields a force per unit of volume, 



■in I 
k ) 



(14a) 



which we express in units of H 2 /8tt^ . Above, the first 
term describes the vortex interaction with the current 
induced by the external field and with the interface, 



F s = — 
' R 



1 



H R 2 



and the second, the vortex-vortex interaction, 



r l - r k 
Ti - r k Y 



(14b) 



(14c) 



The simple Euler method was used to accomplish the 
time integration, but adopting a St value small enough 
to avoid large variations of the vortex positions between 
two consecutive steps. Moreover, the dynamical matrix 
(the Hessian matrix of Gl), whose elements are given by 



d 2 G L 



(15) 



was calculated for the final vortex configuration. In this 
equation, the Greek indexes stand for the components 
of the vector pi, while the Italic indexes are the labels 
for the vortices. The computation of the dynamical ma- 
trix eigenvalues allowed us to tell whether the given state 
was stable or unstable (for a stable state all the dynam- 
ical matrix eigenvalues must be non-negative). Unstable 
states were discarded. 

One difficulty in simulating this system is the fact that 
both Gl and the forces acting on the vortices diverge at 
the disk edge. To overcome this, during the MD sim- 
ulation whenever a vortex was at a distance less than 
£ from the disk edge, it is taken out from the system, 
i.e., this vortex disappears. Therefore, the final num- 
ber of vortices may not be the same as in the beginning. 
This does not lead to any serious concern, since we col- 
lect all the final results from each trial and sort them 
in ascending order of energy. It also allows to compare 
energies of systems containing different number of vor- 
tices for the same external magnetic field and investigate 
which of them correspond to the lower energy, i.e., is the 
ground state. 



III. LOW L -STATES: VORTEX MOLECULES 

In this section we present the results calculated from 
the GL and London theories for low L— states for a thin 
disk of radius R ~ 50£. A comparison between ground 
states in the GL theory and the London approximation 
was done in Ref. I20I for the case of a small disk radius 
(i.e. R — 6£). In that case, it was not possible to study 
multivortex configurations for L— states above L = 14 
since the calculated GL results showed only giant vor- 
tices. Moreover, above L = 26 the disk was driven to the 
normal state. In the present case, multivortex configu- 
rations are obtained for much higher L— states. This en- 
abled us to compare large multivortex configurations cal- 
culated by both the GL theory and the London approx- 
imation, and investigate the transition to the Abrikosov 
lattice. 

For L — 1 to L — 9, the lowest energy configurations 
consist of vortices distributed in regular polygons with 
or 1 vortex in the center of the disk. This means that 
not many meta-stable states are close to the ground state, 
which makes the job of finding low energy configurations 
easier. In the London limit, this reduces Eqs. (|12afl - l|12d[) 
to a simple form, which depends on only one free param- 
eterj2Sy&i i.e., the radius of the ring which circumscribes 
the polygon, /3 r i ng . The minimization problem is then 
straightforward. We also obtained the positions of the 
vortex ring by finding the roots of 
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1 -r 2 



2r 2 
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<H 1 + r 4 



— 2r 2 cos 4> n 



0, (16) 



which follows from the balance of forces acting on each 
vortex [cf. Eq (|14f> ] ^ Here N is the number of vortices 
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(a) 



(b) 





(c) 



(d) 




(e) 



FIG. 1: Vortex configurations for L = (3) and H = 0.007 (a), 
L = (6) and H = 0.01 (b), L = (1, 5) and H = 0.01 (c), and 
L = (1, 6) and Ho = 0.011 (d). The black lines are the contour 
lines of ^(r)! 2 , whereas the white circles indicate the position 
of the vortices according to the London approximation. In (e) 
we show the phase of the order parameter for the L = (6) state 
at Ho ~ 0.022 obtained from the GL equations (on the left) 
and from the London approximation (on the right). 



on the ring (or the number of sides of the polygon) , r = 
Pring/R, 0n = 2nn / N , h = H R 2 /2 and the plus (minus) 
sign should be taken if there is one (zero) vortex in the 
center of the disk. 

A comparison between the calculated GL and London 
vortex config uratio ns is depicted in Fig, ffl The states 
L = 3 (Fig. [LTaH, L = (6) (Fig M , L = (1,5) 
(Fig. 1(c) I, and L = 7 (Fig. 1(d) I were obtained at 



H a = 0.007, H = 0.01, H = 0.01, and H = 0.011, 
respectively. The vortex positions practically coincide 
for the same configurations in both theories. 

The agreement between the vortex positions yielded by 
both theories (at Hq <C H-di) is related to the fact that 
the phase of the order parameter, 9, is well described as 
the imaginary part of the complex function 



- W0) 2 



EL 
R 



Ho 
4 



(R 2 -p 2 ), (17) 



for sufficiently small magnetic fields^ where £ = pe 1 ^ — 
x + iy is the representation of the vector p in the com- 
plex (— plane. But (d/ n 2 )Re{fl} is simply the stream- 
line function [cf. Eq. ©] calculated in the London limit. 
That is greatly responsible for the fact that p r ing is vir- 
tually the same in both theories for Ho <C H C 2 . Fig. 1 (e) 
presents the numerically calculate phase of the order pa- 
rameter (left) and the theoretical one obtained from the 
imaginary part of Eq. I|17fl (right) for the state with L = 6 
at H = 0.022. 

The dependence of /3 r ing upon Hq is shown in Fig. ffib) 
obtained within the GL (squares) and the London limit 
(solid line) for the L = 1, (2), (3), (4), (5), (6), (1,6), 
(1,7) states. Both theories predict the same values of 
/9 r in g and, thus, the same stable configurations, as func- 
tion of Ho- Fig. E£b) also shows the radial position over 
which a given regular polygon configuration is not stable 
(dashed lines) as function of Ho (obtained in the London 
limit). The magnetic field in which the stable and unsta- 
ble Pring lines start to departure from each other (open 
circles) mark the onset of stability for each configuration. 
The unstable /O r i ng lines merge to 



Ri 1 



HoR 2 ' 



(18) 



This is simply the position after which the attractive 
force acting on each vortex by its own image becomes 
larger than the force produced by the shielding currents 
(which pulls the vortices inside), as can be easily demon- 
strated from Eqs. I|14a|l and (|14b|l for one vortex. It is 
also important to take into account the vortex interac- 
tion with the disk edge for sufficiently low fields. This 
can be noticed from the difference between the stable 
Pring and the dotted lines in Fig. EJb), which depicts the 
position at which the respective regular polygon config- 
uration would sit if there were no vortex images [from 
Eq. (|16f) in the absence of vortex images, p r i ng would be 
given by yj (N ± 1)/Hq, where the + (— ) sign should be 
considered for one (zero) vortex in the center] . 

The free energies within the GL (thick lines) theory 
and the usual London limit (dashed line) are depicted 
in Fig. Efa) for L = — > 8 as a function of the ap- 
plied magnetic field Hq. The energy calculated within 
the London limit (with a = 1) starts to departure from 
the GL results as soon as L = 1. This is mainly due to 
the fact that the usual London theory neglects the spa- 
tial variation of l^ 2 !- When the magnetic field increases, 
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5 10 - 15 20 25 

FIG. 2: (a): The GL (thick lines) and the improved London 
(thin lines) free energies as function of the applied field Ho for 
low L~ states. The L = (1,5) state has slightly lower energy 
than the L = (6) state, as seen in the inset, where the lines 
and the squares show the difference between the L — (6) and 
L = (1,5) energies in the London limit and in the GL the- 
ory, respectively. The usual London energy (where we added 
— 1) is also depicted (dashed lines) for comparison. The solid 
circles show the points at which the usual London energy pre- 
dicts a transition from L to L + 1. (b): the GL (open squares) 
and London (solid lines) radial position of the vortices in the 
ring (pring) as function of the magnetic field for the L = 1, 
(2), (3), (4), (5), (6), (1,6), (1,7) states. The arrow indicates 
the direction of increasing L. For each configuration (in the 
London limit) the vortex position at which the vortex ring is 
unstable (dashed lines) and the onset field from which stabil- 
ity occurs (open circles) are depicted. The radial positions of 
the vortex ring when the boundary induced 'vortex images' 
are neglected are shown by the dotted lines for comparison. 



the ground state changes by the addition of one vortex, 
i.e, L = 0^1^2...^8 (for the London limit these 
transitions are marked by the filled circles). For disks 
with small radius the GL theory predicts that L = 2 — > 6 
states do not have a vortex in the center of the diski^2i24 
Such a central vortex appears in the L = 7 — > 9— states. 
In contrast, for the present large disk case (i? = 50£), 
the GL theory and the London approximation yield five 



vortices arranged in a regular pentagon with one in the 
center of the disk for L = 6. The state with six vortices 
in a regular hexagon has a slightly higher energy [the 
difference in energy is depicted in the inset of Fig.|2Ja)]. 

In an effort to remedy the differences in the energy be- 
tween the GL and the usual London results we considered 
the contribution of the vortex cores energies to the Lon- 
don energy. As long as vortices are well separated and 
H <C 1 (1^1 2 ~ 1 far from the vortex cores), Eq. il4"c|) 
can be approximately given by the London energy. In 
this limit the depreciation of" | 1 2 around the vortex cores 
can be approximated by the superposition of some func- 
tion which varies from to 1 within \p — pj| ~ £. Such 
extensions of the London theory were previously consid- 
ere d 49 ' 50 for infinite superconducting systems, e.g., by 

using |*| 2 = \p — pi\ 2 I (\p — pi\ 2 + 2^ close to the 

core of the vortex at pi. We used this expression into 
Eq. I|4bjl in the limit that vortices are far apart, i.e., for 
low L values, where we can make use of the superposition 
principle. First, Eq. (|4b|) can be written as 
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d 2 p.(19) 



Close to the cores, 1 



*| 2 

3/2 



2/ [\p-Pi\ +2) and 



V|*| = 2/ Hp - Pi\ + 2) (remembering that £ = 1 

in our units). Since these expressions rapidly approach 
zero, we approximated the integration over the disk area 
in Eq (|19|) by the sum of integrations around of the vortex 
cores. This yields 



-1 + L 



R 2 ' 



(20) 



We added the above value of <7 C ore to the London energy, 
Ql, assuming that the vortex core have a radius V2£, 
which yields a = \[2 in g coro . The resulting improved 
London energies are presented in Fig. GJa) by thin lines 
for the L = 1, (2), (3), (4), (5), (1, 5), (1,6), and (1,7) 
states. The agreement between this improved London 
theory with the GL results is very good. Such extension 
of the London limit yields the region over which each con- 
figuration is the ground state with much more confidence 
than the usual London limit. 

In the above approximation for Q COTe the depreciation 
of the order parameter near the disk edge was neglected. 
In order to have an estimate of the behavior of | 'J | 2 close 
to p = R, we may consider the first GL equation written 
as 



V 2 * 



1*1 



n 2 



o. 



(21) 



with boundary conditions _ R = and p ■ H\ p R — 

0. Notice that II = V6 A = v — A automatically 
satisfies its boundary condition if u is considered within 
the London limit [cf. Eq. (JHJ)]. For a giant vortex state, 
I* | 2 is radially symmetric, and u = <j>L/ p. For a regular 
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polygon vortex configuration and after averaging v along 
the angular direction, one finds v — <f)LQ(p — /Q r jng)/Pj 
where 8(x) is the Heaviside step function. Therefore, 
one may approximate the superconducting electron den- 
sity by l^appl 2 f» 1 — (L/p — pH /2) 2 inside a ring with 
internal radius, R\, taken somewhat larger than p r ing 
and external radius smaller than R — £ (since the term 
V 2x f in the first GL equation becomes more important 
within distances of £ close to the disk edge). |* a pp| 2 
is minimal at p = R and consequently we can use its 
value at the boundary in order to estimate when the de- 
preciation of |^| 2 close to the edge becomes important 
(notice that the actual |^| 2 is higher close to the disk 
edge than our approximate result, since there is a cor- 
rection of order V 2 \f r /\I r , with V 2 ^ > 0, in this region). 
We found that a 5% depreciation in |*]/ a pp(i?)| 2 (which 
would mean |*(i?)| 2 > 0.95), requires that H w 0.009 
for L = 0, Ho ~ 0.0098 for L = 1, H « 0.0106 
for L = 2, H Q « 0.0114 for L = 3, H « 0.0122 for 
L = 4, H w 0.013 for L = 5, H w 0.0138 for L = 6, 
i?o = 0.0146 for L = 7 and ff = 0.0154 for L = 8, 
which are magnetic field values well above the respec- 
tive regions where each of theses states are the ground 
state. Also the order parameter depreciation close to the 
disk edge results in a less rapid increment of the energy 
of each L— state compared with the energy found within 
the London limit. But for Ho <C H C 2, such difference 
only becomes pronounced at fields well above the mag- 
netic field region over which the respective L— state is the 
ground state. Nevertheless, the depreciation of the order 
parameter close to the edges is important if one wishes 
to understand the entry and exit of vortices in a finite 
system. 



IV. HIGH L— STATES: ABRIKOSOV LATTICE 

For large values of the vorticity an Abrikosov lattice 
appears in the interior of the disk. In this section we will 
consider Hq > 0.03 and investigate from which value of 
L the Abrikosov lattice start to occupy a substantial area 
in the center of the disk. 

One difficulty which arises when studying the high 
L— states is due to the fact that the energy difference be- 
tween two different L— states and the energy difference 
between distinct configurations with the same L can be 
comparable and very small. This is illustrated in Fig. 
where the energy of the meta-stable states obtained in the 
London limit at Ho = 0.1 and at Ho = 0.2 are shown. 
For instance, the difference between the two lowest en- 
ergy L = 110 and L — 112 states is less than 10~ 4 . At 
Hq = 0.1 (Ho = 0.2) we found that a vortex configuration 
with L — 111 (L = 234) has the lowest London energy. 
Of course it is always possible that configurations with 
lower London energies have not been reached by our sim- 
ulations (due to the fact that we have a finite number of 
trials, i.e., we made typically 1000 trials). Nevertheless, 
the small difference in the energies give us confidence that 
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FIG. 3: Energies of the meta-stable states (L = 109 — * 115 
and L = 226 — > 237) obtained from simulations within the 
London limit at H /H c2 = 0.1 (left) and H /H c2 = 0.2 
(right). The energy difference between two different L— states 
is comparable to the energy difference between distinct con- 
figurations at the same L— state. 



some of these configurations are at least very close to the 
true ground state within the London limit. Moreover, at 
such high L values, it is expected that the energy yielded 
by the London approximation differs considerably from 
the more realistic results obtained from the GL theory. 

In order to circumvent the limitations of the London 
limit in the calculation of the energy, meta-state states 
are also investigated within the GL theory. In this case, 
the correct contribution to the energy from the spatial 
dependence of ^(r)! 2 is taken into account. Again, the 
question concerning whether the calculated configura- 
tions are the true ground states can be addressed, since 
it is possible that the numerical solution of Eqs. 11121) 
becomes trapped in some local minimum. Nonetheless, 
thermal fluctuations are always present in experiments, 
making some excited states close to the ground state 
available for the system. In addition, there is the al- 
ready mentioned fact that the difference between ener- 
gies in these high L— states is very small. Therefore, the 
achievement of the ground state is not crucial for the 
present study. 

Although the London limit fails to give the precise 
value of the vortex system energy at high L, we expect 
that the vortex positions obtained within such approach 
are in good accordance with the GL results (cf. SectionlTTl 
and Ref.HJ, at least at fields up to H ~ 0.23^i There- 
fore the stability of the 'London' configurations within 
the framework of the GL theory was investigated by 
solving Eqs. and (J2J starting from the given Lon- 
don configuration (usually the ones with lowest energy). 
By using this procedure, we found that the L ~ 110 and 
L ~ 230 configurations, as obtained within the London 
theory, are also stable within the GL formalism. The cal- 
culated GL energies of such configurations are very close 
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FIG. 4: Superconducting electron density for L = 44, 56, 
64, 79, 88, 104, 109, and 229 obtained at H = 0.04, 0.05, 
0.06, 0.07, 0.08, 0.09, 0.10 and 0.20, respectively. The white 
lines depict the Delaunay triangulation for the vortex core 
positions. 



to other GL configurations with the same vorticity, the 
relative difference lying typically between 10~ 4 — 10~ 5 . 
Such values are usually 5 to 10 times smaller than the 
relative energy difference between the L and L + l lowest 
energy states. 

Some of the stable configurations at H = 0.04, 0.05, 
0.06, 0.07, 0.08, 0.09, 0.10 and 0.20 are depicted in Fig.H 
for L = 44, 56, 64, 79, 88, 104, 109, and 229, respectively. 
From the Delaunay triangulation performed for the core 
positions, it can be seen that a triangular vortex config- 



uration in the center of the disk starts to appear as L 
increases. First, for L = 64 and L = 79, an hexagonal 
vortex arrangement starts appearing in the center of the 
disk. Such arrangement begins occupying a larger area 
with increasing vorticity. For L > 100 the Abrikosov lat- 
tice is already present in a considerable region inside the 
disk. 

For the high L— states there is a competition between 
the ring-like structure imposed by the disk geometry, and 
the hexagonal lattice favored by the vortex-vortex inter- 
action. As a result, rings are generally formed close to 
the disk edge while an Abrikosov lattice is present in the 
center of the disk. In order to study the configurations 
obtained within the GL theory, we computed the posi- 
tions of the vortex cores from the calculated ^(f )| ■ 

First we investigate the ringlike structure near the disk 
edge by computing the number of vortices, A, and the av- 
erage density of vortices, < a{p) >— N(p)/2-KpAp, as a 
function of p. These quantities can suggest where ringlike 
structures are formed, since N(p) (as well as < a(p) >) 
should present sharp peaks where ringlike patterns exist. 
For this purpose we divided the disk radius into radial 
strips of length Ap — 1.25£ and counted the number of 
vortices in each of these pieces. N(p) and < a(p) > are 
shown in Fig. for L = 109, L = 229, L = 473 and 
L = 717 at H = 0.1, H = 0.2, H = 0.4 and H = 0.6, 
respectively. The L — 109 and L — 229 were obtained by 
solving the GL equations starting with the L = 110 and 
L = 230 less energetic configurations calculated within 
the London limit. We also plotted the respective configu- 
rations inside each figure. To help the visualization, rings 
were drawn for the two outermost shells and a Delaunay 
triangulation was made for the vortices in the interior of 
these rings. Clearly, both N(p) and < a(p) > have one 
sharp peak near the disk edge, an indication of a ring-like 
structure. This can be observed in the vortex configura- 
tions since the outermost vortices are almost perfectly 
aligned in a ring. For the L = 109 state, both N(p) 
and < u(p) > have additional peaks in the interior of 
the disk. As the vortex configuration also indicates, this 
could be interpreted as a second (deformed) outer ring 
with a somewhat deformed hexagonal lattice in the cen- 
ter. For L — 229, it is clear that vortices are distributed 
in ring-like structures for the two outermost rings with an 
inner Abrikosov lattice. Similar features are present in 
the other L ~ 110 and L ~ 230 vortex states, i.e., sharp 
peaks near the disk edge are also present in N(p) and 
< a(p) >, indicating two outermost ringlike vortex distri- 
bution with an Abrikosov lattice in the center (again this 
Abrikosov lattice is much better defined for L ~ 230) . It 
is also worth to mention that the two outer peaks present 
at L ~ 110 and L ~ 230 are situated around the same 
values of p for configurations calculated within both the 
GL and the London theories. For example, for L = 109 
the peaks are at p w 35 and p w 43, with an empty region 
around p sa 39 and another for p > 45. Moreover the re- 
gions comprised by the peaks in < a(p) > at p sa 35 and 
p w 43 contain 28 and 33 vortices, respectively. In the 




FIG. 5: Number of vortices N(p) (circles) and the average vortex density < a(p) > (solid line) for L = 109, L = 229, L = 473 
and L = 717 at, respectively, Ho = 0.1, Ho = 0.2, Ho = 0.4 and Ho = 0.6. The respective configurations are depicted in the 
insets. The well-defined peaks close to R = 50£ is indicative of a ringlike structure close to the edge. This is also indicated 
by the configurations in the insets, where we plotted rings for the two outermost shells and the Delaunay triangulation for the 
inner vortices. 



case L = 229 one sharp peak occurs around p ss 46. The 
radial region close to this peak contains 48 vortices, with 
no vortices for p > 47. The radial region around the peak 
at p « 40 has 44 vortices, with the region between these 
two maxima, around p ~ 43, also vortex free. A more 
complete description of the number of vortices in the two 
outer rings is presented in Table [I] Taking the number 
of vortices in the first and second outermost rings for the 
configurations given in this Table, as well as other config- 
urations not shown here with the same vorticity, we find 
that the number of vortices in these shells are around, 
respectively, 33 - 34 and 28 ± 1 for L ~ 110 (50 ± 2 and 
45 ±1 for L ~ 230). 

In Fig. the states L = 473, at H = 0.4 and L = 717, 
at Hq = 0.6, are also depicted. As expected, the peaks 
become broader deep inside the disk, suggesting that the 
ring-like structure smears out as one approaches the cen- 
ter of the disk. In addition, as the value of L increases 
the average density becomes more uniform, but preserv- 
ing at least two sharp peaks near the edge. For L = 473 
and L = 717 the most external ring is situated at p ~ 47 



TABLE I: Number of vortices (N) and approximate radial 
position of the two outer shells (< p >), and the bond-angular 
order factor Go for configurations with lower energy. Here (i) 
means all vortices, except the ones belonging to the outermost 
shell; (ii), vortices not at the two outer rings and (Hi) vortices 
at p < 25. 
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FIG. 6: The density-density correlation function (right) and the probability p(8) to find two adjacent nearest neighbors of a 
given vortex within an angle 9 (left) for L = 109 at H = 0.1, L = 229 at H = 0.2, L = 473 at H = 0.4, and L = 717 at 
Ho ~ 0.6. The dashed, solid, and thin solid lines represent p(9) calculated for vortices (i) not in the outermost ring, (ii) not in 
the two outer rings, and (Hi) at p < 25, respectively. 



and contains 70 and 92 vortices, respectively. Notice that 
the two outer rings have a very different number of vor- 
tices which is quite distinct from the situation of classi- 
cal charges confined by a parabolic potential 4 - where for 
large number of charges the outer rings contain the same 
number of particles, The present situation is between a 
hard wall 5? and a parabolic confinement case. 

We calculated the density-density correlation function 
for the vortices sitting inside the two outermost rings in 
order to help characterize whether a Abrikosov lattice is 
formed away from the disk edge. This quantity is de- 
picted at the right side of Fig. The density-density 
correlation function indicates an hexagonal pattern for 
all these high L— states. Such pattern is well defined 
for L — 109 at Hq — 0.1 and becomes very well de- 
fined for L — 229 at Ho = 0.2. Other configurations 
with L ~ 110 have also an hexagonal pattern as the one 
for L = 109 (but not as sharp). The density-density 
correlation function computed for various configurations 
with L ~ 230 also resembles the one depicted here for 
L = 229. For L = 473 and L = 717 the hexagonal pat- 
tern is also observed, but not as sharp as the one for 
L = 229. Particularly for the L = 717 configuration, the 
density-density correlation function suggests that each 
vortex (inside the two outermost rings) still have coor- 
dination number equals to six, although the hexagonal 
structure considering the farther neighbors is not well de- 
fined. Therefore these two configurations may still have 
local, but not orientational order beyond some few neigh- 
bors. We shall come back to this point later in SectionlVl 



when discussing the defects in the vortex lattice. 

From the density-density correlation function it is also 
worth to compute the typical inter- vortex distance, a v , 
for the vortices forming the Abrikosov lattice. We thus 
obtained a v ~ 8 for L = 109 at Hq = 0.1, a v ~ 5.8 for 
L = 229 at H = 0.2, a v « 4.1 for L = 473 at H = 0.4, 
and a v w 3.4 for L = 717 at H = 0.6. 

In order to better describe how close the system is to 
an Abrikosov lattice we computed the probability distri- 
bution, p(9), to find two adjacent nearest neighbors of 
a given vortex making an angle 9. This probability was 
calculated for three different cases: (i) for all vortices, ex- 
cept the ones at the outermost ring; (ii) for the vortices 
not in the two outer rings, and (Hi) for those vortices 
at p < 25. These probabilities are shown on the left of 
Figs. El We found that p(9) (for all the cases (i) — * (Hi)) 
is maximum close to 60°, which is characteristic of an 
hexagonal lattice. The width of the distribution rapidly 
decreases as L increases from ~ 110 to ~ 230, but in- 
creases as L is further incremented. To be more precise, 
p(9) for the L ~ 110 (not only the L = 109 state which 
is shown) state obtained within the London limit has a 
maximum at 57°, but with < 9 >= 60° for the cases 
(i) — ► (Hi) (<> means average over the vortices included 
in each case, (i), (ii) or (Hi)). The probability distri- 
butions for cases (i) — * (Hi) are not sharp, presenting 
width of about 12° at half of the distribution maximum. 
Other states with L ~ 110 and comparable energy also 
show similar behavior. Such features can be understood 
as the result of the contribution to the p(9) distribution 
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from vortices in the border of the Abrikosov lattice re- 
gion. Since not so many vortices are present in this region 
for L ~ 110, vortices in its border will contribute more 
strongly to the p(9) distribution than for higher L states. 
Such vortices have to adjust themselves to the ring-like 
structure more than the inner vortices and, so, it is likely 
that a few of them may have nearest neighbors within an- 
gles less than 60° or, even, coordination number different 
to six. For L > 200, p(9) is sharply peaked at 9 = 60°, in 
conformity with the density-density correlation function, 
signaling an Abrikosov lattice in the interior of the disk. 

For completeness we also calculate the bond-angular 
order factoa^^, 

Ge = ( TT £ exp(z A, ib 0„) \ , (22) 

\ iV,lb n=l I 

where N n b = 6 is the number of nearest neighbors of a 
given vortex, 9 n is the angle between two segments join- 
ing the given vortex with two adjacent nearest neighbors 
and <> is again the average over the vortices in cases 
(i), (ii) or (Hi). It is clear from Eq. (1221 that Gq = 1 
for an ideal Abrikosov lattice. In table \£\Gq is depicted 
for some of the configurations we obtained (typically the 
configurations with lowest energy). The values found for 
Gq are larger than 0.9 at the region (Hi) for L ~ 230, 
which indicates a configuration very close to an hexag- 
onal lattice. The L ~ 110 states obtained at Hq = 0.1 
have lower Gq, which corroborates our previous analysis 
suggesting that an Abrikosov lattice is formed but not yet 
occupying a large area inside the disk. Again, for L = 473 
and L — 717 Gq is not as large as the one calculated at 
L ~ 230, but is still close or larger than 0.9 in region 
(iii), which indicates that a local orientational hexagonal 
order is present. In fact for such large L— values Gq no 
longer increases and the peak in p(9) is slightly broad- 
ened due to the appearance of grain boundaries in the 
Abrikosov lattice as will be shown in the next section. 



V. HIGH STATES: DEFECTS IN THE 
VORTEX LATTICE 

As a result of the competition between the geometry 
induced ring-like structure near the disk border and the 
hexagonal structure in the center, topological defects in 
the lattice appear in between these two regions (a fea- 
ture also observed in confined classical systems^S). In 
order to study the distribution of these deffects in the 
disk, we applied the Voronoi construction. In an in- 
finite system both the GL theory and the London ap- 
proach predict a coordination number equal to six and 
the Voronoi construction would yield hexagonal unit cells 
for each vortex. In the disk the situation is different, 
vortices near the edge have to adjust themselves to the 
boundary. Therefore, topological defects in the vortex 
lattice will be present. We shall use the term (wedge) 
disclination for vortices which have a closed unit cell in 



the Voronoi construction with coordination number dif- 
ferent from six. This difference is called the topological 
charge of the disclination. Notice that some vortices at 
the outermost shell have open unit cells in the Voronoi 
construction. For such vortices the expected number of 
nearest neighbors should be four. So in order to de- 
fine the topological defects also for these vortices, the 
topological charge there is defined as the number of first 
neighbors minus 4. By such convention it can be shown 
from Euler's theorem^ that the net topological charge 
in a disk equals —6. In addition, dislocations (a bounded 
pair of one + and one — disclinations) may also appear, 
whose net topological charge is null, in order to adjust 
the vortex system to a configuration with lower energy. 

Fig.0shows the Voronoi construction for the L = 109 
(H = 0.1), L = 111 (H = 0.1), L = 234 (H = 0.2), 
L = 229 (H = 0.2), L = 473 (H = 0.4), and L = 717 
(Ho = 0.6). In all of them it is quite clear that an 
Abrikosov vortex lattice is formed inside the disk, as in- 
dicated in previous section, but with the formation of 
topological defects in the configurations. The net topo- 
logical charge for all configurations obtained (including 
the ones not shown here) is always —6, in accordance with 
the Euler theorem^ However the total absolute charge 
can be much larger than 6. Negatively charged disclina- 
tions (vortices with coordination number < 6) are always 
present. Vortices with coordination number > 6 (posi- 
tive topological charge) appear accompanied by negative 
topological charges, leading to the formation of disloca- 
tions. The defects in the vortex configurations are more 
suitable to sit in the disk edge or in the region delimiting 
the Abrikosov lattice and the ring-like structure. Never- 
theless, as L increases, dislocations proliferate and form 
grain boundaries in the region where the hexagonal lat- 
tice appears. This is also the reason why the L = 473 and 
L = 717 states have smaller Gq values and less sharper 
peaks in the p(9) distribution than the lower L states, 
for instance L ~ 230. Such feature is also observed in 
simulations performed by Reefman and Brom2£ consid- 
ering 2000 vortices (although they considered vortices in 
the London limit without interaction with the disk edge) 
and in classical systems of charged particles interacting 
with each other via the Coulomb potential and confined 
to a parabolic potential^ 

Koulakov and ShklovskiiS described the presence of 
dislocations in configurations of classical charged parti- 
cles confined by a parabolic potential as due to two main 
reasons: the inhomogeneity in the density of particles 
and the presence of disclinations. The latter (which is 
always present in an hexagonal arrangement confined to 
a disk) causes a large deformation in the particle con- 
figurations. Dislocations thus appear in order to reduce 
such deformations, eventually decreasing the energy of 
the system. Such effect, also called screening, was previ- 
ously described by Nelson and HalperinS^ when studying 
the melting driven by dislocations in two dimensional 
systems, and is linked to the lack of translational long- 
range order in two-dimensional solid systems (although 
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orientational order is still present)^ These dislocations 
are arranged close to or at the disk edge. The former 
reason induces dislocations in the interior of the disk. In 
Ref. |5f|, it was found that there exists a threshold num- 
ber of particles (which in their case is approximately 700) 
below which dislocations are due mainly from screening 
and, above which, such defects appear due to the inho- 
mogeneity of the particle density. At least qualitative 
similarities exists between such system of charged parti- 
cles and our vortex configurations. Therefore, it is rea- 
sonable to speculate that the same mechanisms which 
drive the appearance of dislocations is also present here. 
Just like in the system of charged particles, dislocations 
are mostly distributed close to and at the disk edge for 
L < 230 and start proliferating in the Abrikosov lattice 
for larger L. 

Finally, in order to further investigate the relation be- 
tween defects in the vortex configurations and the energy 
of the system, we computed the total number of defects 
(the number of the + and — topological charges) in each 
stable configuration obtained within the London frame- 
work. The results are shown in Figs. |$1 for L = 110, 
111 and 112 at H = 0.1 (left) and L = 230, 232 and 
234 at Hq = 0.2 (right). The absolute value of the net 
topological charge is depicted as a solid horizontal line 
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FIG. 8: Number of defects (solid points) versus the energy 
for some of the configurations obtained from the London ap- 
proach. The straight horizontal line is the absolute value of 
the net topological charge. 



and is always equal to six as required by the Euler the- 
orem. The total number of defects - which is directly 
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FIG. 9: Superconducting electron density for H — 0.6 (left) 
and H = 1.02 (right). White to black runs from low to high 
values of |*| 2 . 




FIG. 10: Contour plots of log |»I>| 2 at the center (left) and 
close to the edge of the disk for H — 1. 

related to the number of dislocations in the configura- 
tions - is depicted as points connected by lines. One can 
notice that the total number of defects is not a mono- 
tonic function of the London energy of the configuration. 
Instead, it highly fluctuates. For example, a configura- 
tion free of dislocations (in which only six disclinations 
occur) almost always has a higher energy than, e.g., one 
with a total number of 16 topological charges. This hap- 
pens, for example, for L — 111 at Hq — 0.1 where such 
a configuration with only six disclinations (and no dislo- 
cations) has Gl — 0.1302066, which is 0.5% higher than 
the energy of the lowest energy state, Gl — 0.12958384 
(the Voronoi contruction of the latter configuration is the 
L = 111 state depicted in Fig.0). This indicates that the 
presence, as well as the distribution, of dislocations in the 
vortex configurations plays an important role in lowering 
the energy of such configurations. 

VI. SURFACE SUPERCONDUCTIVITY 

When the external magnetic field approaches Hq = 1 
(or Hq = H C 2 in not normalized units) the vorticity, L, 
becomes large. Inside a thin layer close to the disk edge 



the superconducting electron density, | \I/ 1 2 , is larger than 
in the interior of the diskiS Such a behavior may be un- 
derstood as a result of the superposition of the supercon- 
ducting electron density depreciation close to each vortex 
inside the disk, which is less strong for vortices at the sur- 
face. This already takes place for Hq = 0.6 with L = 717, 
but is highly pronounced at Hq > 1.0. At Hq = 0.6, a 
multivortex state (as was shown in previous figures and 
also on the left of Fig. EJ) is enclosed by this superconduct- 
ing sheath. Within this sheath |^| 2 w 0.75, opposed to 
a maximum of \^\ 2 rs 0.5 between two adjacent vortices. 
Nevertheless, according to the criterion adopted to char- 
acterize the existence of a giant vortex state (l^l 2 < 10~ 4 
in the region between vortices)^ a giant vortex state ap- 
pears at Hq = 1.02. In this state \^>\ 2 < 10~ 4 , except at 
R - 2£ < p < R where 0.2 < |*| 2 < 0.45 (cf. Fig. at 
right). At Hq = 1 the maximum value of |\&| 2 is ~ 10~ 2 
in the region between two adjacent vortex cores, while 
l^l 2 ~ 0.55 at the disk edge. Such a configuration is not 
yet a giant vortex state, although the multivortex state 
in this case is extremely 'dilute'. Possibly H = 1 is close 
to the field in which a giant vortex state decays into a 
multivortex stated Moreover, at this magnetic field the 
depreciation of | , F| 2 close to the vortex cores is different 
whether a vortex sits in the outermost ring or in the inte- 
rior of the disk. This feature is depicted in Fig. 1101 where 
a contour plot of the logarithm of the superconducting 
electron density is shown in the center of the disk (at left) 
and close to the edge (at right). 



VII. CONCLUSIONS 

We investigated the magnetic field dependence of vor- 
tex states in thin disks with large radius. The nonlinear 
GL equations, as well the London approximation were 
used to obtain stable vortex configurations. Although 
both methods lead, for small fields, to similar vortex con- 
figurations, the energies are different. This is the reason 
for the failure of the London limit to yield the correct 
ground state configuration. For low values of the vortic- 
ity we improved the London approximation by including 
the spatial variation of |\1/| 2 close to the vortex cores, 
which resulted in energies which were very close to those 
of the GL approach. 

Multivortex states were obtained for fields up to Ho w 
H C 2, above which a giant vortex state appears. We in- 
vestigated how the configuration of this multivortex state 
changes as function of the magnetic field. At low mag- 
netic fields (Hq <C 0.1i? C 2) we find vortex configurations 
having ringlike distribution, as expected from symme- 
try considerations. However as the number of vortices 
increases, the vortex-vortex repulsion starts playing a 
larger role and we observed the appearance of an hexag- 
onal lattice. The ringlike structure is replaced by an 
Abrikosov lattice in the center of the disk as soon as the 
field is close to 0.1if C 2, when L ~ 100, but is preserved 
near the edges. For fields larger than 0.1 this Abrikosov 
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lattice becomes even more pronounced compared to the 
ringlike structure. 

The topological defects in the vortex configurations 
and their distribution were also studied. We observed 
two types of defects: (wedge) disclinations and disloca- 
tions. The net topological charge is always —6, as re- 
quired for an hexagonal structure confined to a circular 
geometry. Similar to classical particles confined in radi- 
ally symmetric potentials, we find that these topological 
defects appear mostly close to the edge for L < 230, in 
order to adjust the ringlike structure to the Abrikosov 
lattice. We attribute the presence of dislocations in that 
region due to the screening of disclinations. As L in- 
creases further dislocations start to be spread in the cen- 
ter of the disk and form grain boundaries. 

Surface superconductivity was observed at fields 
around and above 0.6H C 2- This surface superconductiv- 
ity becomes more pronounced as the vorticity increases, 
which resulted in a larger overlap between the vortices. 
We also noticed that the transition from a multivortex to 



a giant vortex state takes place at magnetic fields slightly 
above H c i- Just below the formation of the giant vor- 
tex state, the superconducting electron density presents 
markedly distinct spatial dependence close to the disk 
edge - where the vortex structure starts to coalesce - 
compared to what is observed in the center of the disk. 
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